Dynamics of Viscoplastic Deformation in Amorphous Solids 
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We propose a dynamical theory of low-temperature shear deformation in amorphous solids. Our 
analysis is based on molecular-dynamics simulations of a two-dimensional, two-component noncrys- 
talline system. These numerical simulations reveal behavior typical of metallic glasses and other 
viscoplastic materials, specifically, reversible elastic deformation at small applied stresses, irreversible 
plastic deformation at larger stresses, a stress threshold above which unbounded plastic flow occurs, 
and a strong dependence of the state of the system on the history of past deformations. Micro- 
scopic observations suggest that a dynamically complete description of the macroscopic state of this 
deforming body requires specifying, in addition to stress and strain, certain average features of a 
population of two-state shear transformation zones. Our introduction of these new state variables 
into the constitutive equations for this system is an extension of earlier models of creep in metallic 
glasses. In the treatment presented here, we specialize to temperatures far below the glass transition, 
and postulate that irreversible motions are governed by local entropic fluctuations in the volumes 
of the transformation zones. In most respects, our theory is in good quantitative agreement with 
the rich variety of phenomena seen in the simulations. 
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I. INTRODUCTION 



This paper is a preliminary report on a molecular- 
dynamics investigation of viscoplastic deformation in a 
non-crystalline solid. It is preliminary in the sense that 
we have completed only the initial stages of our planned 
simulation project. The results, however, have led us to 
a theoretical interpretation that we believe is potentially 
useful as a guide for further investigations along these 
lines. In what follows, we describe both the simulations 
and the theory. 

Our original motivation for this project was an inter- 
est in the physics of deformations near the tips of rapidly 
advancing cracks, where materials are subject to very 
large stresses and experience very high strain rates. Un- 
derstanding the dissipative dynamics which occur in the 
vicinity of the crack tip is necessary to construct a sat- 
isfactory theory of dynamic fracture. |J Indeed, we be- 
lieve that the problem of dynamic fracture cannot be 
separated from the problem of understanding the condi- 
tions under which a solid behaves in a brittle or ductile 
manner. ||-g] To undertake such a project we eventually 
shall need sharper definitions of the terms "brittle" and 
"ductile" than are presently available; but we leave such 
questions to future investigations while we focus on the 
specifics of deformation in the absence of a crack. 

We have chosen to study amorphous materials because 
the best experiments on dynamic instabilities in fracture 
have been carried out in silica glasses and polymers. J7]|| 
We know that amorphous materials exhibit both brittle 
and ductile behavior, often in ways that, on a macro- 
scopic level, look very similar to deformation in crystals. 
|9| More generally, we are looking for fundamental prin- 
ciples that might point us toward theories of deformation 
and failure in broad classes of macroscopically isotropic 
solids where thinking of deformation in terms of the dy- 
namics of individual dislocations |g|| is either suspect, 



due to the absence of underlying crystalline order, or sim- 
ply intractable, due to the extreme complexity of such 
an undertaking. In this way we hope that the ideas pre- 
sented here will be generalizable perhaps to some poly- 
crystalline materials or even single crystals with large 
numbers of randomly distributed dislocations. 

We describe our numerical experiments in Section 
Q. Our working material is a two-dimensional, two- 
component, noncrystalline solid in which the molecules 
interact via Lennard- Jones forces. We purposely main- 
tain our system at a temperature very far below the glass 
transition. In the experiments, we subject this material 
to various sequences of pure shear stresses, during which 
we measure the mechanical response. The simulations re- 
veal a rich variety of behaviors typical of metallic glasses 
pTJ| |l3[ and other viscoplastic solids, 14 specifically: re- 
versible elastic deformation at small applied stresses, irre- 
versible plastic deformation at somewhat larger stresses, 
a stress threshold above which unbounded plastic flow 
occurs, and a strong dependence of the state of the sys- 
tem on the history of past deformations. In addition, the 
molecular-dynamics method permits us to see what each 
molecule is doing at all times; thus, we can identify the 
places where irreversible molecular rearrangements are 
occurring. 

Our microscopic observations suggest that a dynam- 
ically complete description of the macroscopic state of 
this deforming body requires specifying, in addition to 
stress and strain, certain average features of a popula- 
tion of what we shall call "shear transformation zones." 
These zones are small regions, perhaps consisting of only 
five or ten molecules, in special configurations that are 
particularly susceptible to inelastic rearrangements in re- 
sponse to shear stresses. We argue that the constitutive 
relations for a system of this kind must include equa- 
tions of motion for the density and internal states of 
these zones; that is, we must add new time-dependent 



state variables to the dynamical description of this sys- 
tem. |15| , |l6| Our picture of shear transformation zones is 
based on earlier versions of the same idea due to Argon, 
Spaepen and others who described creep in metallic alloys 
in terms of activated transitions in intrinsically heteroge- 
neous materials. [|17|-f22fl These theories, in turn, drew on 
previous free- volume formulations of the glass transition 
by Turnbull, Cohen and others in relating the transition 
rates to local free- volume fluctuations. ]20| , |2"3] -p5| None of 
those theories, however, were meant to describe the low- 
temperature behavior seen here, especially the different 
kinds of irreversible deformations that occur below and 
above a stress threshold, and the history dependence of 
the response of the system to applied loads. 

We present theory of the dynamics of shear transfor- 
mation zones in Section pTJ| . This theory contains four 
crucial features that are not, so far as we know, in any 
previous analysis: First, once a zone has transformed 
and relieved a certain amount of shear stress, it cannot 
transform again in the same direction. Thus, the system 
saturates and, in the language of granular materials, it 
becomes "jammed." Second, zones can be created and 
destroyed at rates proportional to the rate of irreversible 
plastic work being done on the system. This is the ingre- 
dient that produces a threshold for plastic flow; the sys- 
tem can become "unjammed" when new zones are being 
created as fast as existing zones are being transformed. 
Third, the attempt frequency is tied to the noise in the 
system, which is driven by the strain rate. The stochastic 
nature of these fluctuations is assumed to arise from ran- 
dom motions associated with the disorder in the system. 
And, fourth, the transition rates are strongly sensitive 
to the applied stress. It is this sensitivity that produces 
memory effects. 

The resulting theory accounts for many of the features 
of the deformation dynamics seen in our simulations. 
However, it is a mean field theory which fails to take 
into account any spatial correlations induced by interac- 
tions between zones, and therefore it cannot explain all 
aspects of the behavior that we observe. In particular, 
the mean- field nature of our theory precludes, at least for 
the moment, any analysis of strain localization or shear 
banding. 



II. MOLECULAR-DYNAMICS EXPERIMENTS 

A. Algorithm 

Our numerical simulations have been performed in the 
spirit of previous investigations of deformation in amor- 
phous solids p6|-p9[. We have examined the response 
to an applied shear of a noncrystalline, two-dimensional, 
two-component solid composed of either 10,000 or 20,000 
molecules interacting via Lennard- Jones forces. Our 
molecular dynamics (MD) algorithm is derived from a 
standard "NPT" dynamics scheme M, i.e. a pressure- 



temperature ensemble, with a Nose-Hoover thermostat, 
pl|i33[ and a Parinello-Rahman barostat f3J,|35|] modi- 
fied to allow imposition of an arbitrary two-dimensional 
stress tensor. The system obeys periodic boundary con- 
ditions, and both the thermostat and barostat act uni- 
formly throughout the sample. 

Our equations of motion are the following: 
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Here, r„ and p n are the position and momentum of 
the n'th molecule, and F„ is the force exerted on that 
molecule by its neighbors via the Lennard- Jones inter- 
actions. The quantities in brackets, e.g. [i] or [a], are 
two-dimensional tensors. T is the temperature of the 
thermal reservoir; V is the volume of the system (in this 
case, the area), and N is the number of molecules. Tki n is 
the average kinetic energy per molecule divided by Boltz- 
mann's constant Isb- [c] is the externally applied stress, 
and [a av ] is the average stress throughout the system 
computed to be 
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where F^ m is the i'th component of the force between 
particles n and m; r 3 nm is the j'th component of the vec- 
tor displacement between those particles; and V is the 
volume of the system. L is the locus of points which de- 
scribe the boundary of the simulation cell. While (2.5) 
is not directly relevant to the dynamics of the particles, 
keeping track of the boundary is necessary in order to 
properly calculate intermolecular distances in the peri- 
odic cell. 



The additional dynamical degrees of freedom in (2.1 



2.5|) are a viscosity £, which couples the system to the 
thermal reservoir, and a strain rate, [e] via which the 
externally applied stress is transmitted to the system. 
Note that [i] induces an affine transformation about a 
reference point Ro which, without loss of generality, we 
choose to be the origin of our coordinate system. In a 
conventional formulation, [a] would be equal to —P[I], 
where P is the pressure and [I] is the unit tensor. In that 
case, these equations of motion are known to produce the 
same time-averaged equations of state as an equilibrium 
NPT ensemble. [B0l By instead controlling the tensor [a], 
including its off-diagonal terms, it is possible to apply a 
shear stress to the system without creating any preferred 
surfaces which might enhance system-size effects and in- 
terfere with observations of bulk properties. The applied 
stress and the strain-rate tensor are constrained to be 



symmetric in order to avoid physically uninteresting ro- 
tations of the cell. Except where otherwise noted, all of 
our numerical experiments are carried out at constant 
temperature, with P = 0, and with the sample loaded 
in uniform, pure shear. 

We have chosen the artificial time constants tt and Tp 
to represent physical aspects of the system. As suggested 
by Nose [pTl] , Tp is the time for a sound wave to travel an 
interatomic distance and, as suggested by Anderson [f36[ , 
rp is the time for sound to travel the size of the system. 



B. Model Solid 

The special two-component system that we have cho- 
sen to study here has been the subject of other investiga- 
tions |37H3y] primarily because it has a quasi-crystalline 
ground state. The important point for our purposes, 
however, is that this system can be quenched easily into 
an apparently stable glassy state. Whether this is ac- 
tually a thermodynamically stable glass phase is of no 
special interest here. We care only that the noncrys- 
talline state has a lifetime that is very much longer than 
the duration of our experiments. 

Our system consists of molecules of two different sizes, 
which we call "small" (S) and "large" (L). The interac- 
tions between these molecules are standard 6-12 Lennard- 
Jones potentials: 



U a p(r) = 4e Q/3 
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where the subscripts a, (3 denote S or L. We choose the 
zero-energy interatomic distances, a a /3, to be 



a S s = 2sin(— ), a LL = 2sin(-), 
10 5 



o S £ = i; (2. 



with bond strengths: 



ssl = 1, ess = e-LL 



(2.9) 



For computational efficien cy, we impose a finite-range 
cutoff on the potentials in (2/7) by setting them equal to 
zero for separation distances r greater than 2.5asL- The 



masses are all taken to be equal. The ratio of the num- 
ber of large molecules to the number of small molecules 
is half the golden mean: 



Nl 

N s 



1 + V5 



(2.10) 



In the resulting system, it is energetically favorable for 
ten small molecules to surround one large molecule, or for 
five large molecules to surround one small molecule. The 
highly frustrated nature of this system avoids problems 
of local crystallization that often occur in two dimensions 
where the nucleation of single component crystalline re- 
gions is difficult to avoid. As shown by Langon et al [ B7L 
this system goes through something like a glass transition 
upon cooling from its liquid state. The glass transition 
temperature is 0.3To where kp To = esh- All the simula- 
tions reported here have been carried out at a tempera- 
ture T — 0.001 To, that is, at 0.3% of the glass transition 
temperature. Thus, all of the phenomena to be discussed 
here take place at a temperature very much lower than 
the energies associated with the molecular interactions. 

In order to start with a densely packed material, we 
have created our experimental systems by equilibrating 
a random distribution of particles under high pressure 
at the low temperature mentioned above. After allowing 
the system to relax at high pressure, we have reduced the 
pressure to zero and again allowed the sample to relax. 
Our molecular dynamics procedure permits us to relax 
the system only for times of order nanoseconds, which 
are not long enough for the material to experience any 
significant amount of annealing, especially at such a low 
temperature. 

We have performed numerical experiments on two dif- 
ferent samples, containing 10,000 and 20,000 molecules 
respectively. All of the simulation results shown are from 
the larger of the two samples; the smaller sample has been 
used primarily to check the reliablility of our procedures. 
We have created each of these samples only once; thus 
each experiment using either of them starts with precisely 
the same set of molecules in precisely the same positions. 
As will become clear, there are both advantages and un- 
certainties associated with this procedure. On the one 
hand, we have a very carefully controlled starting point 
for each experiment. On the other hand, we do not know 
how sensitive the mechanical properties of our system 





Molecules 


Shear Modulus 


Bulk Modulus 


2D Poisson Ratio 


Young's Modulus 


Sample 1 


10,000 


9.9 


31 


0.51 


30 


Sample 2 


20,000 


16 


58 


0.57 


50 



TABLE I. Sample Sizes and Elastic Constants 



might be to details of the preparation process, nor do 
we know whether to expect significant sample-to-sample 
variations in the molecular configurations. To illustrate 
these uncertainties, we show the elastic constants of the 
samples in Table 1. The moduli are expressed there in 
units of esL/asL- (Note that the Poisson ratio for a two- 
dimensional system has an upper bound of 1 rather than 
0.5 as in the three-dimensional case.) The appreciable 
differences between the moduli of supposedly identical 
materials tell us that we must be very careful in drawing 
detailed conclusions from these preliminary results. 



C. Simulation Results: Macroscopic Observations 
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In all of our numerical experiments, we have tried sim- 
ply to mimic conventional laboratory measurements of 
viscoplastic properties of real materials. The first of these 
is a measurement of stress at constant strain rate. As 
we shall see, this supposedly simplest of the experiments 
is especially interesting and problematic for us because 
it necessarily probes time-dependent behavior near the 
plastic yield stress. 

Our results, for two different strain rates, are shown in 
Figure [fl. The strain rates are expressed in units propor- 
tional to the frequency of oscillation about the minimum 
in the Lennard- Jones potential, specifically, in units of 
luq = (esL/i f na 2 SL )^ , where m is the particle mass. (The 
actual frequency for the SL potential, in cycles per sec- 
ond, is (3 • 2 1 / 3 /7r) uj = 1.2 o>o-) As usual, the sample 
has been kept at constant temperature and at pressure 
P = 0. At low strain, the material behaves in a linearly 
elastic manner. As the strain increases, the response 
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FIG. 1. Shear stress vs. strain for strainrates of 1CP 4 
(solid lines) and 2 x 1CP 4 (dotted lines). The thicker 
lines which denote the simulation results exhibit both lin- 
ear elastic behavior at low strain and non-linear response 
leading to yield at approximately a s — 0.35. The thin- 
ner curves are predictions of the theory for the two strain 
rates. Strainrate is measured in units of (esL/ma 2 SL )^ . 



FIG. 2. Shear strain (open symbols) vs. time for sev- 
eral applied shear stresses (solid symbols). The stresses 
have been ramped up at a constant rate until reaching a 
maximum value and then have been held constant. The 
strain and stress axes are related by twice the shear modu- 
lus so that, for linear elastic response, the open and closed 
symbols would be coincident. For low stresses the sample 
responds in an almost entirely elastic manner. For inter- 
mediate stresses the sample undergoes some plastic defor- 
mation prior to jamming. In the case where the stress is 
brought above the yield stress, the sample deforms indef- 
initely. Time is measured in units of (ma% L /esL)^ ■ 



becomes nonlinear, and the material begins to deform 
plastically. Plastic yielding, that is, the onset of plastic 
flow, occurs when the strain reaches approximately 0.7%. 
Note that the stress does not rise smoothly and mono- 
tonically in these experiments. We presume that most of 
this irregularity would average out in larger systems. As 
we shall see, however, there may also be more interesting 
dynamical effects at work here. 

In all of the other experiments to be reported here, we 
have controlled the stress on the sample and measured 
the strain. In the first of these, shown in Figure |[ we 
have increased the stress to various different values and 
then held it constant. 

In each of these experimental runs, the stress starts at 
zero and increases at the same constant rate until the de- 
sired final stress is reached. The graphs show both this 
applied stress (solid symbols) and the resulting strain 
(open symbols), as functions of time, for three different 
cases. Time is measured in the same molecular- vibration 
units used in the previous experiments, i.e. in units of 
(ma 2 SL /esL)^ ■ The stresses and strain axes are related 
by twice the shear modulus so that, if the response is 
linearly elastic, the two curves lie on on top of one an- 
other. In the case labelled (A), the final stress is small, 
and the response is nearly elastic. For cases (o) and (□), 
the sample deforms plastically until it reaches some final 
strain, at which it ceases to undergo further deformation 
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FIG. 3. Stress and strain vs. time for one particular 
loading where the stress has been ramped up to a s — 0.25, 
held for a time, and then released. Note that, in addi- 
tion to the shear response, the material undergoes a small 
amount of dilation. 
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FIG. 4. Elastic and inelastic strain vs time for the 
same simulation as that shown in Figure H. The inelastic 
strain is found by subtracting the linearly elastic strain 
from the total strain. Note the partial recovery of the in- 
elastic portion of the strain which occurs during and after 
unloading. 

on observable timescales. (We cannot rule out the possi- 
bility of slow creep at much longer times.) In case (O), 
for which the final stress is the largest of the three cases 
shown, the sample continues to deform plastically at con- 
stant stress throughout the duration of the experiment. 
We conclude from these and a number of similar experi- 
mental runs that there exists a well defined critical stress 
for this material, below which it reaches a limit of plastic 
deformation, that is, it "jams," and above which it flows 
plastically. Because the stress is ramped up quickly, we 
can see in curves (□) and (O) of Figure || that there is a 
separation of time scales between the elastic and plastic 



responses. The elastic response is instantaneous, while 
the plastic response develops over a few hundred molec- 
ular vibrational periods. To see the distinction between 
these behaviors more clearly, we have performed experi- 
ments in which we load the system to a fixed, subcritical 
stress, hold it there, and then unload it by ramping the 
stress back down to zero. In Figure ||, we show this stress 
and the resulting total shear strain, as functions of time, 
for one of those experiments. If we define the the elas- 
tic strain to be the stress divided by twice the previously 
measured, as-quenched, shear modulus, then we can com- 
pute the inelastic strain by subtracting the elastic from 
the total. The result is shown in Figure 0. Note that 
most, but not quite all, of the inelastic strain consists 
of nonrecoverable plastic deformation that persists after 
unloading to zero stress. Note also, as shown in Fig- 
ure ||, that the system undergoes a small dilation during 
this process, and that this dilation appears to have both 
elastic and inelastic components. 

Using the simple prescription outlined above, we have 
measured the final inelastic shear strain as a function of 
shear stress. That is, we have measured the shear strain 
once the system has ceased to deform as in the subcritical 
cases in Figure ||, and then subtracted the elastic part. 
The results are shown in Fig. H. As expected, we see only 
very small amounts of inelastic strain at low stress. As 
the stress approaches the yield stress, the inelastic strain 
appears to diverge approximately logarithmically. 

The final test that we have performed is to cycle the 
system through loading, reloading, and reverse-loading. 
As shown in Figure q, the sample is first loaded on the 
curve from a to b. The initial response is linearly elas- 
tic, but, eventually, deviation from linearity occurs as 
the material begins to deform inelastically. From b to 
c, the stress is constant and the sample continues to de- 
form inelastically until reaching a final strain at c. Upon 
unloading, from c to d, the system does not behave in 
a purely elastic manner but, rather, recovers some por- 
tion of the strain anelastically. While held at zero stress, 
the sample continues to undergo anelastic strain recovery 
from d to e. 

When the sample is then reloaded from e to f, it un- 
dergoes much less inelastic deformation than during the 
initial loading. From f to g the sample again deforms 
inelastically, but by an amount only slightly more than 
the previously recovered strain, returning approximately 
to point c. Upon unloading again from g to h to i, less 
strain is recovered than in the previous unloading from c 
through e. 

It is during reverse loading from i to k that it becomes 
apparent that the deformation history has rendered the 
amorphous sample highly anisotropic in its response to 
further applied shear. The inelastic strain from i to k 
is much greater than that from e to g, demonstrating 
a very significant Bauschinger effect. The plastic defor- 
mation in the initial direction apparently has biased the 
sample in such a way as to inhibit further inelastic yield 




FIG. 5. Final inelastic strain vs. applied stress for 
stresses below yield. The simulation data (squares) have 
been obtained by running the simulations until all defor- 
mation apparently had stopped. The comparison to the 
theory (line) was obtained by numerically integrating the 
equations of motion for a period of 800 time units, the 
duration of the longest simulation runs. 
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in the same direction, but there is no such inhibition in 
the reverse direction. The material, therefore, must in 
some way have microstructurally encoded, i.e. partially 
"memorized," its loading history. 



D. Simulation Results: Microscopic Observations 

Our numerical methods allow us to examine what is 
happening at the molecular level during these deforma- 
tions. To do this systematically, we need to identify 
where irreversible plastic rearrangements are occurring. 
More precisely, we must identify places where the molec- 
ular displacements are non-affine, that is, where they de- 
viate substantially from displacements which can be de- 
scribed by a linear strain field. Our mathematical tech- 
nique for identifying regions of non-affine displacement 
has been described by one of us (MLF) in an earlier pub- 
lication, jlij For completeness, we repeat it here. 

We start with a set of molecular positions and sub- 
sequent displacements, and compute the closest possible 
approximation to a local strain tensor in the neighbor- 
hood of any particular molecule. To define that neigh- 
borhood, we define a sampling radius, which we choose 
to be the interaction range, 2.5asL- The local strain is 
then determined by minimizing the mean square differ- 
ence between the the actual displacements of the neigh- 
boring molecules relative to the central one, and the rel- 
ative displacements that they would have if they were in 
a region of uniform strain e^. That is, we define 

£ 2 (t,At)=£]T[<(i)-r«(t)- 

n i 

£(*« + eij)(ri(t - At) - 4(t - At))]\ (2.11) 



where the indices i and j denote spatial coordinates, and 
the index n runs over the molecules within the interac- 
tion range of the reference molecule, n — being the 
reference molecule. r l n {t) is the i'th component of the 
position of the n'th molecule at time t. We then find the 
£ij which minimizes D 2 by calculating: 



FIG. 6. Stress-strain trajectory for a molecular- 
dynamics experiment in which the sample has been 
loaded, unloaded, reloaded, unloaded again, and then re- 
verse-loaded, all at stresses below the yield stress. The 
smaller graph above shows the history of applied shear 
stress with letters indicating identical times in the two 
graphs. The dashed line in the main graph is the theo- 
retical prediction for the same sequence of stresses. Note 
that a small amount of inelastic strain recovery occurs af- 
ter the first unloading in the simulation, but that no such 
behavior occurs in the theory. Thus, the theoretical curve 
from c though h unloads, reloads and unloads again all 
along the same line. 
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The minimum value of D 2 (t, At) is then the local de- 
viation from affinc deformation during the time interval 
[t — At, t). We shall refer to this quantity as D 2 
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FIG. 7. Intensity plots of -D^ini the deviation from 
affine deformation, for various intervals during two simu- 
lations. Figures (a), (b) and (c) show deformation during 
one simulation in which the stress has been ramped up 
quickly to a value less than the yield stress and then held 
constant. Figure (a) shows deformations over the first 
10 time units, and figure (b) over the first 30 time units. 
Figure (c) shows the same state as in (b), but with D^ iri 
computed only for deformations that took place during 
the preceeding 1 time unit. In Figure (d), the initial sys- 
tem and the time interval (10 units) are the same as in 
(a), but the stress has been applied in the opposite direc- 
tion. The gray scale in these figures has been selected 
so that the darkest spots identify molecules for which 

|-Dmin| ~ 0.5 asL- 
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We have found that -D^ in is an excellent diagnostic for 
identifying local irreversible shear transformations. Fig- 
ure u\ contains four different intensity plots of -D^ in for a 
particular system as it is undergoing plastic deformation. 
The stress has been ramped up to |o\,| = 0.12 in the time 
interval [0,12] and then held constant in an experiment 
analogous to that shown in Figure @. Figure R(a) shows 
Dmin f° r * = 10j At = 10. It demonstrates that the non- 
affine deformations occur as isolated small events. In (b) 
we observe the same simulation, but for t = 30, At = 30; 
that is, we are looking at a later time, but again we con- 
sider rearrangements relative to the inital configuration. 
Now it appears that the regions of rearrangement have a 
larger scale structure. The pattern seen here looks like 
an incipient shear band. However, in (c), where t = 30, 
At = 1, we again consider this later time but look only 
at rearrangements that have occurred in the preceeding 
short time interval. The events shown in this figure are 
small, demonstrating that the pattern shown in (b) is, 
in fact, an aggregation of many local events. Lastly, in 
(d), we show an experiment similar in all respects to (a) 
except that the sign of the stress has been reversed. As 
in (a), t = 10, At = 10, and again we observe small 
isolated events. However, these events occur in different 
locations, implying a direction dependence of the local 
transformation mechanism. 

Next we look at these processes in yet more detail. 
Figure His a close-up of the molecular configurations in 
the lower left-hand part of the largest dark cluster seen 
in Figure uKc), shown just before and just after a shear 
transfornration. During this event, the cluster of one 
large and three small molecules has compressed along 
the top-left/bottom-right axis and extended along the 
bottom-lcft/top-right axis. This deformation is consis- 
tent with the orientation of the applied shear, which is in 
the direction shown by the arrows on the outside of the 
figure. Note that this rearrangement takes place without 
significantly affecting the relative positions of molecules 
in the immediate environment of the transforming region. 
This is the type of rearrangement that Spaepen identifies 
as a "flow defect." [EOJ As mentioned in the introduction, 
we shall call these regions "shear transformation zones." 



III. THEORETICAL INTERPRETATION OF THE 
MOLECULAR-DYNAMICS EXPERIMENTS 



A. Basic Hypotheses 



FIG. 8. Close-up picture of a shear transformation zone 
before and after undergoing transfornration. Molecules af- 
ter transformation are shaded according to their values of 
D m in using the same gray scale as in Figure |jj. The di- 
rection of the externally applied shear stress is shown by 
the arrows. The ovals are included solely as guides for the 
eye. 



We turn now to our attempts to develop a theoretical 
interpretation of the phenomena seen in the simulations. 
We shall not insist that our theory reproduce every detail 
of these results. In fact, the simulations are not yet com- 
plete enough to tell us whether some of our observations 
are truly general properties of the model or are artifacts 
of the ways in which we have prepared the system and 
carried out the numerical experiments. Our strategy will 



be, first, to specify what we believe to be the basic frame- 
work of a theory, and then to determine which specific 
assumptions within this framework are consistent with 
the numerical experiments. 

There are several features of our numerical experi- 
ments that we shall assume are fundamentally correct 
and which, therefore, must be outcomes of our theory. 
These are: (1) At a sufficiently small, fixed load, i.e. un- 
der a constant shear stress less than some value that we 
identify as a yield stress, the system undergoes a finite 
plastic deformation. The amount of this deformation di- 
verges as the loading stress approaches the yield stress. 
(2) At loading stresses above the yield stress, the system 
flows visco-plastically. (3) The response of the system to 
loading is history-dependent. If it is loaded, unloaded, 
and then reloaded to the same stress, it behaves almost 
elastically during the reloading, i.e. it does not undergo 
additional plastic deformation. On the other hand, if it 
is loaded, unloaded, and then reloaded with a stress of 
the opposite sign, it deforms substantially in the opposite 
direction. 

Our theory consists of a set of rate equations describ- 
ing plastic deformation. These include an equation for 
the inelastic strain rate as a function of the stress plus 
other variables that describe the internal state of the sys- 
tem. We also postulate equations of motion for these 
state variables. Deformation theories of this type are in 
the spirit of investigations by E. Hart |lj| who, to the 
best of our knowledge, was the first to argue in a math- 
ematically systematic way that any satisfactory theory 
of plasticity must include dynamical state variables, be- 
yond just stress and strain. A similar point of view has 
been stressed by Rice, pi Our analysis is also influenced 
by the use of state variables in theories of friction pro- 
posed recently by Ruina, Dieterich, Carlson, and others. 

Our picture of what is happening at the molecular level 
in these systems is an extension of the ideas of Turnbull, 
Cohen, Argon, Spaepen and others. |17 H24P3H 25[ These 
authors postulated that deformation in amorphous mate- 
rials occurs at special sites where the molecules are able 
to rearrange themselves in response to applied stresses. 
As described in the preceding chapter, we do see such 
sites in our simulations, and shall use these shear trans- 
formation zones as the basis for our analysis. However, 
we must be careful to state as precisely as possible our 
definition of these zones, because we shall use them in 
ways that were not considered by the previous authors. 

One of the most fundamental differences between pre- 
vious work and ours is the fact that our system is ef- 
fectively at zero temperature. When it is in mechanical 
equilibrium, no changes occur in its internal state be- 
cause there is no thermal noise to drive such changes. 
Thus the shear transformation zones can undergo tran- 
sitions only when the system is in motion. Because the 
system is strongly disordered, the forces induced by large- 
scale motions at the position of any individual molecule 
may be noisy. These fluctuating forces may even look as 



if they have a thermal component. J47| The thermody- 
namic analogy (thermal activation of shear transforma- 
tions with temperature being some function of the shear 
rate) may be an alternative to (or an equivalent of) the 
theory to be discussed here. However, it is beyond the 
scope of the present investigation. 

Our next hypothesis is that shear transformation zones 
are geometrically identifiable regions in an amorphous 
solid. That is, we assume that we could — at least in 
principle — look at a picture of any one of the computer- 
generated states of our system and identify small regions 
that are particularly susceptible to inelastic rearrange- 
ment. As suggested by Fig. || these zones might con- 
sist of groups of four or more relatively loosely bound 
molecules surrounded by more rigid "cages." But that 
specific picture is not necessary. The main idea is that 
some such irregularities are locked in on time scales that 
are very much longer than molecular collision times. 
That is not to say that these zones are permanent fea- 
tures of the system on experimental time scales. On the 
contrary, the tendency of these zones to appear and dis- 
appear during plastic deformation will be an essential 
ingredient of our theory. 

We suppose further that these shear transformation 
zones are two-state systems. That is, in the absence of 
any deformation of the cage of molecules that surrounds 
them, they are equally stable in either of two configu- 
rations. Very roughly speaking, the molecular arrange- 
ments in these two configurations are elongated along 
one or the other of two perpendicular directions which, 
shortly, we shall take to be coincident with the principal 
axes of the applied shear stress. The transition between 
one such state and the other constitutes an elementary in- 
crement of shear strain. Note that bistability is the natu- 
ral assumption here. More than two states of comparable 
stability might be possible but would have relatively low 
probability. A crucial feature of these bistable systems 
is that they can transform back and forth between their 
two states but cannot make repeated transformations in 
one direction. Thus there is a natural limit to how much 
shear can take place at one of these zones so long as the 
zone remains intact. 

We now consider an ensemble of shear transformation 
zones and estimate the probability that any one of them 
will undergo a transition at an applied shear stress a s . 
Because the temperatures at which we are working are 
so low that ordinary thermal activation is irrelevant, we 
focus our attention on entropic variations of the local 
free volume. Our basic assumption is that the transition 
probability is proportional to the probability that the 
molecules in a zone have a sufficiently large excess free 
volume, say AV*, in which to rearrange themselves. This 
critical free volume must depend on the magnitude and 
orientation of the elastic deformation of the zone that is 
caused by the externally applied stress, a s . 

At this point, our analysis borrows in its general ap- 
proach, but not in its specifics, from recent developments 
in the theory of granular materials Eq| where the only 



extensive state variable is the volume il. What follows is 
a very simple approximation which, at great loss of gen- 
erality, leads us quickly to the result that we need. The 
free volume, i.e. the volume in excess of close packing 
that the particles have available for motion, is roughly 



tt- Nv = Nv f , 



(3.1) 



where N is the total number of particles, Vf is the av- 
erage free volume per particle, and vq is the volume per 
particle in an ideal state of random dense packing. In the 
dense solids of interest to us here, Vf -C Vo, and there- 
fore vq is approximately the average volume per particle 
even when the system is slightly dilated. The number 
of states available to this system is roughly proportional 
to (vf/h) N , where h is an arbitrary constant with di- 
mensions of volume — the analog of Planck's constant 
in classical statistical mechanics — that plays no role 
other than to provide dimensional consistency. Thus the 
entropy, defined here to be a dimcnsionlcss quantity, is 



S(n,N) = NlnC^-) 



^ TV In 



n-Nvo 

Nh 



(3.2) 



The intensive variable analogous to temperature is x 



(3.3) 



Our activation factor, analogous to the Boltzmann factor 
for thermally activated processes, is therefore 



p -(AV*/x) Qd e -(AV/«/) i 



(3.4) 



A formula like 
earlier literature. 



appears in various places in the 
5| There is an important dif- 
ference between its earlier use and the way in wh ich we 
are using it here. In earlier interpretations, (3.4) is an 
estimate of the probability that any given molecule has a 
large enough free volume near it to be the site at which 
a thermally activated i rrev ersible transition might occur. 
In our interpretation, (3.4) plays more nearly the role of 
the thermal activation factor itself. It tells us something 
about the configurational probability for a zone, not just 
for a single molecule. When multiplied by the density of 
zones and a rate factor, about which we shall have more 
to say shortly, it becomes the transformation rate per 
unit volume. 

Note what is happening here. Our system is extremely 
non-ergodic and, even when it is undergoing appreciable 
strain, does not explore more than a very small part of 
its configuration space. Apart from the molecular rear- 
rangements that take place during plastic deformation, 
the only chance that the system has for coming close 
to any state of equilibrium occurs during the quench by 
which it is formed initially. Because we control only the 
temperature and pressure during that quench, we must 
use entropic considerations to compute the relative prob- 
abilities of various molecular configurations that result 
from it. 



The transitions occurring within shear transformation 
zones are strains, and therefore they must, in principle, 
be described by tensors. For present purposes, however, 
we can make some simplifying assumptions. As described 
in Chapter [Ffl, our molecular-dynamics model is subject 
only to a uniform, pure shear stress of magnitude a s and 
a hydrostatic pressure P (usually zero). Therefore, in the 
principal-axis system of coordinates, the stress tensor is: 



w 



-P <Ts 

a, -P 



(3.5) 



Our assumption is that the shear transformation zones 
are all oriented along the same pair of principal axes, and 
therefore that the strain tensor has the form: 



£s £d 



(3.6) 



where e s and Ed are the shear and dilational strains re- 
spectively. The total shear strain is the sum of elastic 
and inelastic components: 

e s =ef+sT. (3.7) 

By definition, the elastic component is the linear response 
to the stress: 



2/i' 



(3.8) 



where fi is the shear modulus. 

In a more general formulation, we shall have to con- 
sider a distribution of orientations of the shear transfor- 
mation zones. That distribution will not necessarily be 
isotropic when plastic deformations are occurring, and 
very likely the distribution itself will be a dynamical en- 
tity with its own equations of motion. Our present anal- 
ysis, however, is too crude to justify any such level of 
sophistication. 

The last of our main hypotheses is an equation of mo- 
tion for the densities of the shear transformation zones. 
Denote the two states of the shear transformation zones 
by the symbols + and — , and let n± be the number den- 
sities of zones in those states. We then write: 

h± =R T n T - R± n± - d (a s e*") n± + C 2 (cr s £*")• 

(3.9) 

Here, the R± are the rates at which ± states transform 
to =F states. These must be consistent with the transition 
probabilities described in the preceding paragraphs. 

The last two terms in ( p.9[ ) describe the way in which 
the population of shear transformation zones changes as 
the system undergoes plastic deformation. The zones 
can be annihilated and created — as shown by the terms 
with coefficients C\ and C-i respectively — at rates pro- 
portional to the rate <J s e" L at which irreversible work is 
being done on the system. This last assumption is simple 



and plausible, but it is not strictly dictated by the physics 
in any way that we can see. A caveat: In certain circum- 
stances, when the sample does work on its environment, 
o s i™ could be negativ e, in which case the annihilation 
and creation terms in (3.9) could produce results which 
would not be physically plausible. We believe that such 
states in our theory are dynamically accessible only from 
unphysical starting configurations. In related theories, 
however, that may not be the case. 

It is important to recognize that the annihilation and 
creation terms in (3.E) are interaction terms, and that 
they have been introduced here in a mean-field approx- 
imation. That is, we implicitly assume that the rates 
at which shear transformation zones are annihilated and 
created depend only on the rate at which irreversible 
work is being done on the system as a whole, and that 
there is no correlation between the position at which the 
work is being done and the place where the annihilation 
or creation is occurring. This is in fact not the case as 
shown by Fig. M(b), and is possibly the weakest aspect 
of our theory. 

With the preceding definitions, the time rate of change 
of the inelastic shear strain, e* n , has the form: 



V z Ae [R+n + -R_ri-] 



(3.10) 



where V z is the typical volume of a zone and Ae is the 
increment of local shear strain. 



B. Specific Assumptions 

We turn now to the more detailed assumptions and 
analyses that we need in order to develop our general 
hypotheses into a testable theory. 

According to our hypothesis about the probabilities of 
volume fluctuations, we should write the transition rates 



in (3.9) in the form: 



R± = Rq exp 



AV*{±a s 



v f 



(3.11) 



The prefactor Rq is an as-yet unspecified att empt fre- 
quency for these transformations. In writing (3.11), we 
have used the assumed symmetry of the system to note 
that, if AV*((T S ) is the required excess free volume for a 

(-1 > — ) transition, then the appropriate free volume for 

the reverse transition must be AV*(—a s ). We adopt the 
convention that a positive shear stress deforms a zone in 

such a way that it enhances the probability of a (H ► — ) 

transition and decreases the probability of a (— — > +) 
transition. Then AV* (a s ) is a decreasing function of a s . 
Before going any further in specifying the ingredients 
of Rq, AV*, etc., it is useful to recast the equations of 
motion in the following form. Define 



ntot 



n& = U— — jij 



(3-12) 



C(a s 



S(<T S ) = - 



exp 



exp 



AV*(a s ) 

Vf 



AV*(a s ) 



+ exp 



exp 



AV*(-a s ) 

v f 



AV*(-a s ) 



v f 



(3.13) 



(For convenience, and in order to be consistent with later 
assumptions, we have suppressed other possible argu- 
ments of the functions C(a s ) and S(a s ).) Then (3.11) 
becomes: 



Ro V z Ae 



n tot S(tj s ) - n A C(a s ) 



The equations of motion for ua and ntot are: 



«A 



2s< 



V z Ae 



a s n A 



an c 



and 



n to t 



2 cr s i™ i in,,. 
V z Aea 



where a and noo are defined by 

2 



Ci 



V, Ae rioo a 



C 2 



1 



V z Aea 



(3.14) 



(3.15) 



(3.16) 



(3.17) 



oo is the stable equilibrium 
'" remains positive, a is a 



From ( 3.16| ), we see that n 

value of ntot so long as a s i 

characteristic stress that, in certain cases, turns out to 

be the plastic yield stress. As we shall see, we need only 

the above form of the equations of motion to deduce the 

existence of the plastic yield stress and to compute some 

elementary properties of the system. 

The interesting time-dependent behavior of the sys- 
tem, however, depends sensitively on the as-yet unspec- 
ified ingredients of these equations. Consider first the 
rate factor Rq. Our zero-temperature hypothesis implies 
that Rq should be zero whenever the inelastic shear rate 
e| n and the elastic shear rate if = a s /2a both vanish. 
Accordingly, we assume that 



R 2 v 1 ' 2 [(if) 2 + (iTf 



,1/4 



(3.18) 



and 



where v is a constant that we must determine from the 
numerical data. Note that v contains both an attempt 
frequency and a statistical factor associated with the mul- 
tiplicity of trajectories leading from one state to the other 
in an active zone. |49| 

We can offer onl y a speculative justification for the 
right-hand side of ( 3.18| ). The rearrangements that oc- 
cur during irreversible shear transformations are those in 
which molecules deviate from the trajectories that they 
would follow if the system were a continuous medium 
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undergoing affine strain. If we assume that these devia- 
tions are diffusive, and that the affine deformation over 
some time interval scales like the strain rate, then the 
non-affine transformation rate must scale like the square 
root of the affine rate. (Diffusive deviations from smooth 
trajectories have been observed directly in numerical sim- 
ulations of sheared foams, JlTj but o nly in the equivalent 
of our plastic flow regime.) In (3.18), we further assume 
that the elastic and inelastic strain rates are incoherent, 
and thus write the sum of squares within the brackets. 
In wh at fol lows, we shall not be able to test the valid- 
ity of ( |3.18 ) with any precision. Most probably, the only 
properties of importance to us for present purposes are 
the magnitude of Rq and the fact that it vanishes when 
the shear rates vanish. 

Finally, we need to specify the ingredients of AV* and 
Vf. For AV*, we choose the simple form: 



AV*(a s ) = V *exp(-a s /fl) 



(3.19) 



where V Q * is a volume, perhaps of order the average 
molecular volume vq, and p, has the dimensions of a shear 
modulus. The right- hand side of ( |3.1S ) simply reflects 
the fact that the free volume needed for an activated 
transition will decrease if the zone in question is loaded 
with a stress which coincides with the direction of the 
resulting strain. We choose the exponential rather than 
a linear dependence because it makes no sense for the 
incremental free volume V * to be negative, even for very 
large values of the applied stresses. 

Irreversibility enters the theory via a simple switching 
behavior that occurs when the a s ~ dependence of AV* in 



(I 



19) is so strong that it converts a negligably small rate 
at a s — to a large rate at relevant, non-zero values of 
(t s . If this happens, then zones that have switched in one 
direction under the influence of the stress will remain in 
that state when the stress is removed. 

In the formulation presented here, we consider vj to 
be constant. This is certainly an approximation; in fact, 
as seen in Figure [3|, the system dilates during shear de- 
formation. We have experimented with versions of this 
theory in which the dilation plays a controlling role in 
the dynamics via variations in Vf. We shall not discuss 
these versions further because they behaved in ways that 
were qualitatively different from what we observed in our 
simulations. The differences arise from feedback between 
inelastic dilation and flow which occur in these dilational 
models, and apparently not in the simulations. A sim- 
ple comparison of the quantities involved demonstrates 
that the assumption that Vf is approximately constant 
is consistent with our other assumptions. If we assume 
that the increment in free volume at zero stress must be 
of order the volume of a small particle, V * « vq « 0.3, 
and then look ahead and use our best-fi t value for the 
ratio V */vf « 14.0 (see Section [ill D| , Table ||), we 
find Vf ~ 0.02. Since the change in free volume due 
to a dilational strain Ed is Avf = Ed/ p, where p is the 
number density, and Ed < 0.2% for all shear stresses ex- 
cept those very close to yield, it appears that, generally, 



Avf ss EdVo <C Vf. Even when Ed = 1%, the value 
observed in our simulations at yield, the dilational free 
volume is only about the same as the initial free volume 
estimated by this analysis. 

As a final step in examining the underlying structure of 
these equations of motion, we make the following scaling 
transformations: 

2 ±±^ £ , ^^A; 2~ S A; ^ ee E. (3.20) 



Then we find: 



£ = £?(£, A, A); 



A = 2J r (E, A, A)(1-£A); 



A = 2.F(E, A, A)E(l-A); 



where 



J="(E, A, A) = R Q [A5(E) - AC(E)] 



and: 



C(E) = ~ 2 



5(E) = - 



cxp( _]£ e -^ ) + cxp( _:_i e A S) 



v f 

Vf 



exp^e-^-expt-^O 



Vf 

yi 

Vf 



Here, 



& = -■ £ - 2 M ^ Ae n " 
p,' a 



The rate factor in (3.18) can be rewritten 



R = vi E 2 + £ 



where 



a 
v = — - v. 
2/i 



(3.21) 

(3.22) 
(3.23) 

(3.24) 



(3.25) 



(3.26) 



(3.27) 



(3.28) 



C. Special Steady-State Solutions 

Although in general we must use numerical methods to 
solve the fully time dependent equations of motion, we 
can solve them analytically for special cases in which the 
stress E is held constant. Note that none of the results 



presented in this subsection, apart from ( 3.35 ), depend 
on our specific choice of the rate factor Rq. 

There are two specially important steady-state solu- 
tions at constant E. The first of these is a jammed so- 
lution in which £ = 0, that is ,F(E, A, A) vanishes and 
therefore: 
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A = A^ = AT(E); 
C(£) v ' 



where 



T(E) = 1-2 



1 + cxp (2-2-sinh(AE)j 



(3.29) 



(3.30) 



Now suppose that, instead of increasing the stress at a 
finite rate as we have done in our numerical experiments, 
we let it jump discontinuously — from zero, perha ps — 
to its value E at time t = 0. While S is constant, ( 3.22 ) 



and (3.23) can be solved to yield: 



1 - A(i) _ 1 - EA(t) 
1-A(0) ~ 1-EA(0) : 



where A(0) and A(0) denote the initial valu es of A(i) 
and A (t) respectively. Similarly, we can solve (3.21) and 



( 3.22 ) for £{t) in terms of A(i) and obtain a relationship 
between the bias in the population of defects and the 
change in strain, 



^-<°>-4'°(^§) 



(3.32) 



Combining ( 3.2S ), ( 3.31 ) and ( 3.32J ), we can determine 
the change in strain prior to jamming. That is, for E 
sufficiently small that the following limit exists, we can 
compute a final inelastic strain £f. 



£ f = lim £{t) = 

t—*oo 



£{Q) 



2E ln ^ + L 1 - ET(E) ) 



The right-hand side of ( ^33|) , for £{0) = A(0) = 0, 
should be at least a rough approximation for the inelastic 
strain as a function of stress as shown in Figure |g. 

The preceding analysis is our mathematical descrip- 
tion of how the system jams due to the two-state nature 
of the shear transformation zones. Each increment of 
plastic deformation corresponds to the transformation of 
zones aligned favorably with the applied shear stress. As 
the zones transform, the bias in their population — i.e. 
A — grows. Eventually, all of the favorably aligned zones 
that can transform at the given magnitude and direction 
of the stress have undergone their one allowed transfor- 
mation, A has become large enough to cause T in (3.24) 
to vanish, and plastic deformation comes to a halt. 

The second steady-state is a plastically flo wing solu- 
tion in which £^0butA = A = 0. From (|3.22| ) and 
( |3.23 ) we see that this condition requires: 



A =E ; 



A = l. 



(3.34) 



This leads us directly to an equation for the strain-rate 
at constant applied stress, 



£ = v£ z 



5(E) - -C(E) 



(3.35) 



This flowing solution arises f rom the non-linear annihi- 
lation and creation terms in (3.9). In the flowing state, 
stresses are high enough that shear transformation zones 
are continuously created. A balance between the rate of 
zone creation and the rate of transformation determines 
the rate of deformation. 

Examination of (3.22) and ( [3.23| ) reveals that the 
jammed solution ( 3.29|) i s stable for low stresses, while 
the flowing solution^ 3.34 ) is stable for high stresses. The 
crossover between the two solutions occurs when both 



( |3.29 ) and (3.34) are satisfied. This crossover defines the 



(3.31) yield stress E y , which satisfies the condition 



=-=7"(2v). 

Zjy 



(3.36) 



argument of the logarithm in (3.32) 



Note that the 
diverges at E 
(2VJ/v f )siiih(AEy) > 1, ( |3.36| ) implies that E 



Note also that, so long 



as 
1. 



This inequality is easily satisfied for the parameters dis- 
cussed in the following subsection. Thus the dimensional 
yield stress a y is approx imat ed accurately by a in our 
original units defined in (3.20). 



D. Parameters of the Theory 



There are five adjustable system parameters in our the- 



(3.33) ory: a, V z Aen c 



Vg/vf, and y,. In addition, we 



must specify initial conditions for £, A, and A. For 
all cases of interest here, £ (0) = A(0) = 0. However, 
A(0) = n tot (0) / Hoc is an important parameter that char- 
acterizes the as-quenched initial state of the system and 
which remains to be determined. 

To test the validity of this theory, we now must find out 
whether there exists a set of physically reasonable values 
of these parameters for which the theory accounts for 
all (or almost all) of the wide variety of time-dependent 
phenomena seen in the molecular-dynamics experiments. 
Our strategy has been to start with rough guesses based 
on our understanding of what these parameters mean, 
and then to adjust these values by trial and error to fit 
what we believe to be the crucial features of the exper- 
iments. We then have used those values of the parame- 
ters in the equations of motion to check agreement with 
other numerical experiments. In adjusting parameters, 
we have looked for accurate agreement between theory 
and experiment in low-stress situations where we expect 
the concentration of active shear transformation zones to 
be low; and we have allowed larger discrepancies near and 
above the yield stress where we suspect that interactions 
between the zones may invalidate our mean-field approx- 
imation. Our best-fit parameters are shown in Table [Fj]. 



12 



Parameter 


Value 


a 


0.32 


V z Ae jioo 


5.7% 


V 


50.0 


V?/«/ 


14.0 


M 


0.25 


nt o t(0)/noo 


2.0 



TABLE II. Values of parameters for 
simulation data 



comparison to 



The easiest parameter to fit should be a because it 
should be very nearly equal to the yield stress. That is, 
it should be somewhere in the range 0.30-0.35 according 
t o th e data shown in Fig. Q Note that we cannot use 
(3.33) to fit the experimental data near the yield point 
because both the numerical simulations and the theory 
tell us that the system approaches its stationary state 
infinitely slowly there. Moreover, we expect interaction 
effects to be important here. The solid curve in Fig. @ 
is the theoretically predicted strain found by integrating 
the equations of motion for 800 time units, the duration 
of the longest of the simulation runs. The downward ad- 
justment of ex, from its apparent value of about 0.35 to its 
best-fit value of 0.32, has been made on the basis of the 
latter time-dependent calculations plus evidence about 
the effect of this parameter in other parts of the theory. 

Next we consider V z Ae rioo , a dimensionless parameter 
which corresponds to the amount of strain that would oc- 
cur if the density of zones were equal to the equilibrium 
concentration (ntot = ^oo) and if all the zones trans- 
formed in the same direction in unison. Alternatively, if 
the local strain increment Ae is about unity, then this 
parameter is the fraction of the volume of the system 
that is occupied by shear transformation zones. In either 
way of looking at this quantity, our best-fit value of 5.7% 
seems sensible. 

The parameter v is a rate which is roughly the prod- 
uct of an attempt frequency and a statistical factor. The 
only system-dependent quantity with the dimensions of 
inverse time is the molecular vibrational frequency, which 
we have seen is of order unity. Our best-fit value of 50 
seems to imply that the statistical factor is moderately 
large which, in turn, implies that the shear transforma- 
tion zones are fairly complex, multi- molecule structures. 
Lacking any first-principles theory of this rate factor, 
however, we cannot be confident about this observation. 

Our first rough guess for a value of Vq /vf comes from 
the assumption that AV* must be about one molecular 
volume in the absence of an external stress, and that Vf 
is likely to be about a tenth of this. Thus our best-fit 
value of 14.0 is reassuringly close to what we expected. 

The parameter p,, a modulus that characterizes the 
sensitivity of AV* to the applied stress, is especially in- 
teresting. Our best-fit value of 0.25 is almost two orders 
of magnitude smaller than a typical shear modulus for 



these systems. This means that the shear transforma- 
tions are induced by relatively small stresses or, equiv- 
alently, the internal elastic modes within the zones are 
very soft. This conclusion is supported quite robustly 
by our fitting procedure. Alternative assumptions, such 
as control by variations in the average free volume Vf 
discussed earlier, produce qualitatively wrong pictures of 
the time dependent onset of plastic deformation. 

Finally, we consider A(0) = n to t(0) /noc, the ratio of 
the inital zone density to the equilibrium zone density. 
This parameter characterizes the transient behavior as- 
sociated with the initial quench; that is, it determines the 
as-quenched system's first response to an applied stress. 
We can learn something about this parameter by looking 
at later behavior, i.e. the next few segments of a hystere- 
sis loop like that shown in Figure pi. If, as is observed 
there, the loop narrows after the first leg, then we know 
that there was an excess of shear transformation zones 
in the as-quenched system, and that this excess was re- 
duced in the initial deformation. An initial excess means 
A(0) > 1, consistent with our best-fit value of 2.0. 



E. Comparisons between Theory and Simulations 

We now illustrate the degree to which this theory can 
— and cannot — account for the phenomena observed in 
the numerical experiments. 

Figure summarizes one of the principal successes of 
the theory, specifically, its ability to predict the time- 
dependent onset of plastic deformation over a range of 
applied stresses below the yield stress. The solid lines 
in the Figure show the shear strains in three different 
simulations as functions of time. In each simulation the 
stress is ramped up at the same controlled rate, held 
constant for a period of time, and then ramped down, 
again at the same rate. In the lowest curve the stress 
reaches a maximum of 0.1 in our dimensionless stress 
units {s.sl/<^sl)i m the middle curve 0.2, and in the 
highest 0.3. The dashed lines show the predictions of 
the theory. The excellent agreement during and after 
the ramp-up is our most direct evidence for the small 
value of jl quoted above. The detailed shapes of these 
curves at the tops of the ramps, where a s drops abruptly 
to zero, provide some qualitativ e sup port for our choice 
of the rate dependence of Rq in ( 3.18 ). As shown in Fig- 
ure and discussed in the preceding subsection, the final 
inelastic strains in these ramp-up experiments are also 
predicted adequately by the theory. 

The situation is different for the unloading phases of 
these experiments, that is, during and after the periods 
when the stresses are ramped back down to zero. The 
theoretical strain rates shown in Figure |9| vanish abruptly 
at the bottoms of the ramps because our transformation 
rates become negligably small at zero stress. In the two 
experimental curves for the higher stresses, however, the 
strain continues to decrease for a short while after the 
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FIG. 9. Strain vs. time for simulations in which the 
stress has been ramped up at a controlled rate to stresses 
of 0.1, 0.2 and 0.3, held constant, and then ramped down 
to zero (solid lines) . The dashed lines are the correspond- 
ing theoretical predictions. 

stresses have stopped changing. Our theory seems to rule 
out any such recovery of inelastic strain at zero stress; 
thus we cannot account for this phenomenon except to 
remark that it must have something to do with the ini- 
tial state of the as-quenched system. As seen in Figure ||, 
no such recovery occurs when the system is loaded and 
unloaded a second time. 

In Figure g, we compare the stress-strain hysteresis 
loop in the simulation (solid line) with that predicted by 
the theory (dashed line) . Apart from the inelastic strain 
recovery after the first unloading in the simulation, the 
theory and the experiment agree well with one another 
at least through the reverse loading to point k. The 
agreement becomes less good in subsequent cycles of the 
hysteresis loop, possibly because shear bands are forming 
during repeated plastic deformations. 

In the last of the tests of theory to be reported here, we 
have added in Figure [j] two theoretical curves for stresses 
as functions of strain at the two different constant strain 
rates used in the simulations. The agreement between 
theory and experiment is better than we probably should 
expect for situations in which the stresses necessarily rise 
to values at or above the yield stress. Moreover, the va- 
lidity of the comparison is obscured by the large fluctu- 
ations in the data, which we believe to be due primarily 
to small sample size. 

Among the interesting features of the theoretical re- 
sults in Figure [l] are the peaks in the stresses that occur 
just prior to the establishment of steady-states at con- 
stant stresses. These peaks occur because the internal 
degrees of freedom of the system, specifically A(£) and 
A(t), cannot initially equilibrate fast enough to accomo- 
date the rapidly increasing inelastic strain. Thus there 
is a transient stiffening of the material and a momentary 
increase in the stress needed to maintain the constant 



strain rate. This kind of effect may in part be the ex- 
planation for some of the oscillations in the stress seen 
in the experiments. In a more speculative vein, we note 
that this is our first direct hint of the kind of dynamic 
plastic stiffening that is needed in order to transmit high 
stresses to crack tips in brittle fracture. The orders of 
magnitude of the time scales are even roughly the same. 
The strain rates used here, and those that occur near the 
tips of brittle cracks, are both of the order of 10 7 per 
second. 



IV. CONCLUDING REMARKS 

The most striking and robust conclusion to emerge 
from this investigation, in our opinion, is that a wide 
range of realistic, irreversible, viscoplastic phenomena oc- 
cur in an extremely simple molecular-dynamics model 
- a two-dimensional, two-component, Lennard- Jones 
amorphous solid at essentially zero temperature. An al- 
most equally striking conclusion is that a theory based 
on the dynamics of two-state shear transformation zones 
is in substantial agreement with the observed behavior of 
this model. This theory has survived several quantitative 
tests of its applicability. 

We stated in our Introduction that this is a preliminary 
report. Both the numerical simulations and the theoreti- 
cal analysis require careful evaluation and improvements. 
Most importantly, the work so far raises many important 
questions that need to be addressed in future investiga- 
tions. 

The first kind of question pertains to our molecular- 
dynamics simulations: Are they accurate and repeatable? 
We believe that they are good enough for present pur- 
poses, but we recognize that there are potentially impor- 
tant difficulties. The most obvious of these is that our 
simulations have been performed with very small sys- 
tems; thus, size effects may be important. For example, 
the fact that only a few shear transforming regions are 
active at any time may account for abrupt jumps and 
other irregularities sometimes seen in the simulations, 
e.g. in Figure pi. We have performed the simulations 
in a periodic cell to eliminate edge effects. We also have 
tried to compare results from two systems of different 
sizes, although only the results from the larger system 
are presented here. Unfortunately, comparisons between 
any two different initial configurations are difficult be- 
cause of our inability, as yet, to create reproduceable 
glassy starting configurations (a problem which we shall 
discuss next). However, we have seen qualitatively the 
same behavior in both systems, and assume that phe- 
nomena which are common to both systems can be used 
as a guide for theoretical investigations. 

As noted in Section IIB and in Table , our two systems 
had quite different elastic moduli. (Remarkably, their 
yield stresses were nearly identical. It would be interest- 
ing to learn whether this is a repeatable and/or physically 
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important phenomenon.) The discrepancy between the 
elastic properties of the two systems leads us to believe 
that, in future work, we shall have to learn how to con- 
trol the initial configurations more carefully, perhaps by 
annealing the systems after the initial quenches. Unfor- 
tunately, straightforward annealing at temperatures well 
below the glass transition is not yet possible with stan- 
dard molecular-dynamics algorithms, which can simulate 
times only up to about a microsecond for systems of this 
size even with today's fastest computers. Monte Carlo 
techniques or accelerated molecular-dynamics algorithms 
may eventually be useful in this effort. |50|-|5^] An alter- 
native strategy may be simply to look at larger numbers 
of simulations. 

By far the most difficult and interesting questions, 
however, pertain to our theoretical analysis. Although 
Figures and [| provide strong evidence that irreversible 
shear transformations are localized events, we have no 
sharp definition of a "shear transformation zone." So far, 
we have identified these zones only after the fact, that is, 
only by observing where the transformations are taking 
place. Is it possible, at least in principle, to identify zones 
before they become active? 

One ingredient of a better definition of shear transfor- 
mation zones will be a generalization to isotropic amor- 
phous systems in both two and three dimensions. As we 
noted in Chapter III, our functions n±(t) should be ten- 
sor quantities that describe distributions over the ways 
in which the individual zones are aligned with respect 
to the orientation of the applied shear stress. We be- 
lieve that this is a relatively easy generalization; one of 
us (MLF) expects to report on work along these lines in 
a later publication. 

Our more urgent reason for needing a better under- 
standing of shear transformation zones is that, with- 
out such an understanding, we shall not be able to find 
first-principles derivations of several, as-yet purely phe- 
nomenological, ingredients of our theory. It might be 
useful, for example, to be able to start from the molecu- 
lar force constants and calculate the p aram eters Vq and 
p, that occur in the activation factor ( [3.19] ). These pa- 
rameters, however, seem to have clear physical interpre- 
tations; thus we might be satisfied to deduce them from 
experiment. In contrast, the conceptually most challeng- 
ing and important terms are the rate factor in (3.18) and 
the annihilation and creation terms in (3.9), where we do 



not even know what the function al for ms ought to be. 

Calculating the rate factor in ( 3.18 ), or a correct ver- 
sion of that equation, is clearly a very fundamental prob- 
lem in nonequilibrium statistical physics. So far as we 
know, there are no studies in the literature that might 
help us compute the force fluctuations induced at some 
site by externally driven deformations of an amorphous 
material. Nor do we know how to compute a statisti- 
cal prefactor analogous, perhaps, to the entropic factor 
that converts an activation energy to an activation free 
energy. E3 We do know, however, that that entropic fac- 
tor will depend strongly on the size and structure of the 



zone that is undergoing the transformation. 

As emphasized in Chapter III, the annihilation and 
creation terms in ( |3.9| ) describe interaction effects. Even 
within the framework of our mean-field approximation, 
we do not know with any certainty what these terms 
should be. Our assumption that they are proportional to 
the rate of irreversible work is by no means unique. (In- 
deed, we have tried other possibilities in related investi- 
gations and have arrived at qualitatively similar conclu- 
sions.) Without knowing more about the nature of the 
shear transformation zones, it will be difficult to derive 
such interaction terms from first principles. 

A better understanding of these interaction terms is 
especially important because these are the terms that 
will have to be modified when we go beyond the mean- 
field theory to account for correlations between regions 
undergoing plastic deformations. We know from our sim- 
ulations that the active zones cluster even at stresses far 
below the plastic yield stress; and we know that plastic 
yield in real amorphous materials is dominated by shear 
banding. Thus, generalizing the present mean-field the- 
ory to one which takes into account spatial variations 
in the densities of shear transformation zones must be a 
high priority in this research program. 

Finally, we return briefly to the questions which moti- 
vated this investigation: How might the dynamical effects 
described here, which must occur in the vicinity of a crack 
tip, control crack stability and brittle/ductile behavior? 
As we have seen, our theoretical picture of viscoplasticity 
does allow large stresses to be transmitted, at least for 
short times, through plastically deforming materials. It 
should be interesting to see what happens if we incorpo- 
rate this picture into theories of dynamic fracture. 
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